Bending mode fluctuations and structural stability of graphene nanoribbons 
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We analyze the thermal fluctuations of a narrow graphene nanoribbon. Using a continuum, 
membrane- like model we study the height-height correlation functions and the destabilization modes 
corresponding to two different boundaries conditions: ribbons which are fixed or free on the edges. 
For the first situation, the thermal spectrum has a gap and the correlations along the ribbon decay 
exponentially. Thermal fluctuations produce only local perturbations of the flat situation. However, 
the long range crystalline order is not distorted. For free edges, the situation changes as thermal 
excitations are gapless. The low energy spectrum decouples into a bulk and an edge excitation. 
The bulk excitation tends to destabilize the crystalline order producing an homogeneous rippling. 
Furthermore, we associate the edge mode with a precluding perturbation leading to scrolled edges, 
as seen in suspended graphene samples. 

PACS numbers: 61.48.Gh, 63.22.Rc, 81.07.Bc 
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I. INTRODUCTION 

The interplay between lattice deformations and elec- 
tron dynamics will be an important ingredient to under- 
stand and control the electronic properties of graphene 
future devices. On one side, an external strain applied to 
graphene produces a pseudo-magnetic field whose effect 
has been first theoretically predicted^ and then experi- 
mentally determined^. This could be the starting point of 
a field called straintronics, namely the control of the elec- 
tronic properties by applying mechanical strain. On the 
other hand, the intrinsic corrugation observed since the 
early experiments in graphene suspended samples affects 
the electron mobility. Moreover, the fluctuations of this 
corrugation, called flexural phonons, have been proposed 
to be the source of the intrinsic limit in the electronic 
mobility of graphene suspended samples^. Certainly, the 
control of these corrugations is an important question 
to address in order to improve the electronic mobility in 
graphene. 

When the dimensionality is reduced, height fluctua- 
tions are amplified due to the known tendency to in- 
stabilities in low dimensions. We expect that thick 
ribbons with quasi-one-dimensional geometry will have 
stronger thermal fluctuations than two dimensional sys- 
tems. These fluctuations will have important effects on 
the electronic transport and their mechanism should be 
identified in order to control and manage the electronic 
properties of graphene nanoribbons. 

It is the goal of the present paper to study the thermal 
excitations of graphene nanoribbons. We take as a start- 
ing point a continuum model, allowing us to account for 
the long-wavelenght acoustic phonons. Our focus is to 
understand how different boundary conditions affect the 
vibrational modes and how these vibrations perturbate 
the flatness of the static case. We analyze these points 
by calculating the out-of-plane flexural phonons and the 
height-height correlation functions of two different situa- 



tions: clamped and free edges. 

This paper is organized as follows: in Section II we 
introduce the model by taking a continuum limit of a 
lattice bending system. We also discuss how the appro- 
priate boundary conditions can be taken into account. 
In Section III we present a general formalism based in 
a path integral to obtain the correlation functions. In 
Section IV and V we obtain the out-of-plane phononic 
spectrum, the correlation functions, and analyze their 
consequences. Finally, in Section VI we give our conclu- 
sions and perspectives. 



II. THE MODEL AND THE BOUNDARY 
CONDITIONS 



Single- and few-layer graphene are systems of atomic- 
scale thickness. As such, a continuum elastic theory of 
thick plates is not directly applicable. However, mechan- 
ical properties, formation of ripples and the phonon spec- 
trum as the basis of the electron-phonon interaction, are 
well described by the elastic energy form of thick plates. 
The clue to understand this fact is that the bending rigid- 
ity of graphene does not arise from compressions and di- 
latations of the continuum medium bounded by free sur- 
faces. Therefore, the bending rigidity parameter cannot 
be obtained from the elastic parameters of the medium; 
instead, it is an independent quantity^. It is thought that 
the bending rigidity in graphene is due to the bond-angle 
and bond-order term associated with the dihedral angles 
of the underlying C-C interactions^. 

This distinction has special significance in the presence 
of edges, as is the case of the ribbons we are interested on 
in this work. To make the discussion concrete, we start 
from a simplified lattice bending model which has been 
introduced in the studies of membranes^. The model 
Hamiltonian is: 
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where rii is the unit normal vector at the ith site 
of the lattice and j is its nearest neighbor. We use R 
as the bending rigidity parameter in the lattice model. 
For small deviations from the flat configuration, we can 
parametrize the vector position of a point on the surface 
as: 



r{x,y) = (x,y,h(x,y)). 



(2) 



This is called a Mongc representation. h(x, y) is the 
height variable and the normal unit vector becomes: 



n(x,y) = 



e z - Vh 



(-d x h, -d y h, 1) 



(3) 



where e z is the unit vector in the z direction. To obtain 
a continuum limit of Hamiltonian (fTJ) we undertake a 
gradient expansion. We obtain: 
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continuum 



<i,j> 



lim -k J2 (n* - nj) 2 = | / d 2 x |Vn| 



<»j> 



d 2 x [(d 2 x hf + (d 2 y h) 2 + 2(d x d v h) 2 ] , (4) 



where we have neglected terms of order 0(h 2 ) such as 
(V 2 h) 2 in the denominator of Eq. Q and the one in the 
integration measure. Explicit geometric contributions of 
Eq. (j4]) can be seen if we write it down in the following 

way: 



d 2 x [(Vh) 2 -2Det(didjh)] 



(5) 



The first term is proportional to the square of the mean 
curvature and the last to the Gaussian curvature, both 
written in the lowest order of a gradient expansion in 
h(x,y). In terms of these curvatures, Eq. ([5]) is known 
as the Helfrich form of the bending energy of a liquid 
membrane^. The Gaussian curvature is a total derivative 
term and was neglected in previous studies of the stability 
of graphene membranes. However, it plays an important 
role in the search for the appropriate boundary conditions 
of our ribbon geometry. 

If, instead of Eq. (fTJ), we had begun with the elasticity 
problem of a thin plate, the bending energy would have 
assumed the same expression as Eq. ((5]^. However, in 
this last case, the bending rigidity would not have been 
an independent parameter but a function of the clastic 



modulus of the plate as k = x ^_ a -±y being I its thick- 
ness, Y the Young's modulus and a the Poisson ratio. 
Also, the Gaussian curvature of the second term of Eq. 
([5]) would have been multiplied by (1 — a). Regarding the 
discussion at the beginning of the section, we see that 
there is a formal connection between the theory of thin 
clastic plates and the one of a two dimensional mem- 
branes. The last can be obtained from the former by 
taking the bending rigidity as an independent parameter 
and by setting a = 0. 

In addition to the surface term corresponding to the 
Gaussian curvature, there is another total derivative term 
coming from the first part of Eq. ([5]), which can be writ- 
ten as: 
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h (di + dl + 2d 2 x d 2 ) h + 



d x (d x hd 2 x h - hd 3 x h - 2hd 2 d x h) + 
d y (d y hd 2 h - hd 3 y h + 2d x hd v d x h) 



(0) 



Until now we have not specified neither the integra- 
tion domain nor the physically boundary condition for 
our problem. We consider a long and narrow ribbon 
of width W and length L running along the y direc- 
tion. The domain of integration is therefore given by 
-W/2 x < W/2 and -L/2 ^ y < L/2 with W < L. 
Periodic boundary conditions are taken in the y direc- 
tion. Therefore, the surface term corresponding to the 
last line of Eq. ^ vanishes. We finally obtain the fol- 
lowing bending energy for the ribbon: 



Jy / dx[h(di + dt J + 2d 2 x d 2 )h] 
dy\ d x h(x,y)[d 2 x h{x,y)] - 



h(x, y) [d 3 x h(x, y) + 2d 2 d x h(x, y)] (7) 

J x =±W/2 

The last symbol means that the term into braces is the 
difference between this expression evaluated at x = ^7- 
and at x = — -y-. Therefore, the integral runs along the 
edges of the ribbon. A cancellation of this term could 
take place in two different situations. Consider a ribbon 
with its edges clamped along the y direction. In this case 
we have: 



W 

h(x = ± Y ,y) = 
W 

d x h(x = ± — ,y) = 0. 



(8) 



The terms multiplying h(x = ±^-,y) and d x h(x = 
±y,y) can be interpreted as the force and the torque 
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on the edge of the ribbon^. Setting these terms to zero 
means having free edges, and the boundary conditions 
are then: 



{d 3 x + 2d 2 y d x )h(x = ±™,y) = 
d 2 h(x = ±^-,y) = 0. 



(9) 



Note that those are the situations (clamped and free) 
considered in Ref. to obtain the phonon dispersion re- 
lations of a nanoribbon. In this reference, the theory of 
an elastic thin plate is used and the following boundary 
conditions are considered for free edges: 



This situation is quite similar to the one analyzed in 
Ref. to solve the classical dynamics for the out-of-plane 
normal modes. We define the dimcnsionless variables 



q m = Wq m , A = 
lutions of Eq. (12 



XW 2 , x 
can be written as: 



W and V = w- The sa 



with 



Pi(m,n) = ±\ q 2 n ± A. 



(13) 



(14) 



{dl + (2 - <j)d 2 d x ) h(x 



w 



(d 2 x + *d*)h(x = ±—,y) 



(10) 



For a — 0, Eq. ([T0|) is the same as Eq. (|9]) as we al- 
ready remarked. So far we have not discussed about the 
clastic stretching energy of the ribbon. This is because, 
in the harmonic approximation, the in-plane modes that 
correspond to the stretching energy are decoupled from 
the out-of-plane ones, allowing us to study these situa- 
tions separately. In the next section we will compare the 
" height-height" correlation function and the mean square 
of the height corresponding to the two different boundary 
conditions of Eqs. © and ©• 



III. GENERAL FORMALISM FOR 
CORRELATION FUNCTIONS 

For both clamped and free boundary conditions, the 
surface term in Eq.([5]) vanishes and the partition function 
of the system can be written as a path integral of the 
form: 



d.r 



Vhe 



h (d^ + d^ + 2d 2 d 2 )h 



integrating over all the paths that fulfill the boundary 
condition ([8]) or (J9j. It is convenient to expand the path 
in the basis of the eigenfunctions of the operator O. Due 
to the periodic boundary condition in the y direction,we 
can separate its y dependence. The eigenfunctions as- 
sume the form: f% l (x)e'" ImV ', where q m = and 
are the eigenfunctions of the following problem: 



[q 4 m + d A x - 2q 2 m d 2 x ] /" (x) = X 2 /£(*). (12) 



After replacing Eq. (fl3|) in the boundary conditions ([8]) 
or ([9]) we obtain a linear 4x4 problem whose solutions give 
A as a function of q, and the coefficients di to construct 
the normalized eigenfunctions. We will elaborate this 
route in the next section. 

Once we have solved this problem, the expression for 
h(x,y) can be developed as: 



h(x,y) 



f n (f) p iqmV 



(15) 



and with the corresponding change of variables, 
the path integral (fTTj) will run now on the coeffi- 
cients a mn . The height-height correlation function 
(h(xi, yi)h(x2, §2)) can be obtained as usual by adding 
a source term to the exponent of the form a mi „f mi „ 
and then taking the derivative with respect to e m , n - The 
result is: 



(h(xi,yi)h(x 2 ,y 2 )) 



kT _ c i(yi-y 2 )g„ 

kL ^ 

m,n 

/£(S1)/£(S2) 



(16) 



with h(x,y) = b!^£^JLL_ j n the previous equation we 
assumed that the eigenfunctions are normalized in such 
( n ) a way that f \ | | 2 dx = 1. 



IV. CORRELATION FUNCTIONS FOR 
CLAMPED RIBBONS 

Imposing the conditions of Eq. (J5J over h(x,y) ex- 
pressed in terms of f^(x) according to Eq. (fl"3"j) . we ob- 
tain a 4x4 matrix Ai(X, q) multiplying the vector of co- 
efficients (di, c?2, d^, cLi), whose result must be zero. Ex- 
plicitly: 
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where a 



Qm + -Vi 



COS 

—b sin 
—6 sin 



and b 



MS] \ 

.2. 
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sin 
— sin 
& cos 
— & cos 



(17) 



A, 



9 



E, are real 



variables. By requiring that Det[.M] = we obtain the 
values of X m ,n for each q m . The index n = 0,1,2,3... 
used so far enumerates each one of the dispersion curves, 
which are shown in Fig. [TJ 



1250 



1000 - 



750 - 



I a- 



500 



:5o 




FIG. 1. Dispersion curves given by the functions A(g) for 
the clamped ribbon. We show the first seven branches of the 
spectrum which, in fact, has an infinite number of them. In 
the inset we show a zoom of the low energy spectrum for the 
first two branches. 

These boundary conditions allow us to simplify the 
determinant of the matrix M. , which can be factorized in 
two terms: 



Det[.M] = -8 (a cosh[^] sin[^] - b cos[^] sinh[^]) 

(6 cosh[~] sin[^] + a cos[^] sinh[|]). (18) 

Under these circumstances the matrix M. is reduced 
to a simpler form and the eigenfunctions fm(x) can be 
explicitly calculated. When the first term vanishes, we 
obtain the n = 1,3,5... branches of the spectrum with 
the odd eigenfunctions, which read: 



fm( X ) ~ C n 



cos[6/2] 
cosh[a/2] 



cosh[a:r] + cos[6a;] I , (19) 



and when the second term becomes zero we obtain the 
n = 0, 2, 4... ones with the even eigenfunctions, which can 
be written as: 



fm{ x ) — @n 



2 ^ smh[ax] + sin [6a;] ) . (20) 



sinh[a/2] 



The quantities C n , as remarked at the end of sec- 
tion Mil represent normalization constants. The plots 
for (fm( x )) 2 with n = 0, 1, 2 and q m = 6ir are shown in 
Fig. El 



1.5 

1 

0.5 



n=0 branch 
n=l branch 
n=2 branch 




FIG. 2. Square of the normalized eigenfunctions f£i(x) for the 
first three branches of the spectrum in the clamped ribbon. 
The calculations are made for q — 6ir. 

As has been remarked in Ref. |9|, there is a gap in the 
spectrum and the zero energy mode does not exist for 
q m = 0. This is related to the fact that global transla- 
tions are not allowed because the ribbon is gripped by the 
edges. The gap in the first branch behaves as A ~ 
(in the original units) approaching the zero value for the 
infinite square sheet. We expect that height-height cor- 
relations at different points will decay exponentially and 
this is indeed the case. In Fig. [3] we show the value 
of ^(ft,(0.25, y)h(0.25, 0)) running along the y direction 
and evaluated numerically from Eq. (|16p . The contribu- 
tions of the first three branches are shown. As the gap 
increases as long as we go to branches with higher en- 
ergy, the contributions of the corresponding correlations 
become increasingly smaller. 

We see a rapid decay of the correlations in a distance of 
the order of IF. In fact we can estimate the characteris- 
tic correlation length with the following approximations. 
The first branch of Fig. [1] can be fitted well by a function 
of the form Ao(g) — \J ao + a-iq 2 + a2<z 4 with ao = 500, 
ai = 24 and 02 = 0.972. If we neglect the weak de- 
pendency of the eigenfunctions on q rn in Eq. (|16p . the y 
dependency of the correlation is given by the following 
Fourier transform: 



(h(xi,y)h(x 2 ,0)} 



+°° dq 



2tt a + a\q 2 + a 2 q 



(21) 
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FIG. 3. Height-height $ {h(0.25, y)h(0.25, 0)) correlation 
as a function of the distance in the long direction, for the 
clamped ribbon. The contributions of the three first branches 
are shown separately. The dashed line represents the approx- 
imation given by Eq. ()22|) . 



which can be analytically solved, giving: 

(h(x u y)h(x 2 ,0)) s f m { Xl )f m {x 2 )e-^ 

(a sin (q R y) + f3 cos (quy)) , (22) 

where a = 0.00499, (3 = 0.00271 and q R +iqi = 2.273+ 
iA. 185 is a zero of the denominator of Eq. ([2T]) . The decay 
of the correlation is clearly dominated by the exponential 
term. Its characteristic scale i.e. the correlation length, 
is £ = W/ 4.185 (in the original units). 

We see that it is possible to control the extension of 
the height-height correlation by changing the width of 
the ribbon. If we associate this thermal fluctuation with 
the rippling, this results imply that the characteristic size 
of the rippled region grows linearly with the width of the 
ribbons. Let us discuss how the rippling changes as we 
move in the x-direction. We show in Fig. @] the value of 
{h 2 (x,y)). 

The dominant contribution coming from the first 
branch produces a maximum distortion at the center of 
the ribbons. The other branches produce periodic distor- 
tions according to the shape of the eigenfunctions f^(x), 
as shown in Fig. [5J The number of nodes is exactly n + 2 
including those at the edges. 



V. CORRELATION FUNCTIONS AND 
STABILITY OF FREE EDGE RIBBONS 

In a recent paper— one of the authors has studied the 
thermal fluctuations and the stability of a graphene rect- 
angle. Periodic boundary conditions were used in both 
directions in order to avoid edge effects. It has been 
found, from both analytical calculations and Montecarlo 
simulations, that there is a critical relation between the 
width and the length of the system (called B^d-id) in 
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FIG. 4. Mean square of the height ^(h(x,y) 2 } as a function 
on x, the distance to the center, for the clamped ribbon. We 
show the contributions of the first three branches. 

which the dependence of (h 2 (x,y)) with the system size 
changes. When the width decreases and the relation with 
the length becomes smaller than Bid-xd-, thermal oscil- 
lations are more pronounced than in the opposite case. 
The system behaves as one-dimensional and there is a 
higher tendency to a crumpled instability than in square 
samples. 

Moreover, the results of Ref. [l^ do not take into ac- 
count possible excitations of the edges. To account for 
this kind of excitations in a more realistic description, 
we analyze the out-of-plane modes of a ribbon with free 
boundary conditions, as given by Eq. Following sim- 
ilar steps as in the previous section, but with the corre- 
sponding boundary conditions, we obtain the dispersion 
relations of the out-of-plane modes which are shown in 
Fig. [5j Differently than in the clamped ribbon, as trans- 
lational invariance is not broken, there is no gap but two 
acoustic branches with zero energy at q = 0. 
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FIG. 5. Dispersion relations giving the functions A(q) for the 
free ribbon. We show the first seven branches of the spectrum. 
In the inset we show a zoom of the low energy spectrum for 
the two first branches. 
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Qualitatively, wc find a similar spectrum as in Ref.i 
where a finite Poisson ratio is taken into account and 
the boundary conditions given in Eq. (|10j) are used. 
There are, however, some important differences which 
appear with a careful inspection of the two lowest energy 
branches and their corresponding eigenfunctions. First 
of all, note that a function of the form h(x,y) = Ce iqv , 
in which C is a constant, fulfills condition © but not 
(|T0|) . This is a solution of the eigcnproblcm ([12]) with a 
constant eigenfunction and dispersion relation A (q) = q 2 , 
and this is precisely the eigenfunction corresponding 
to the lowest energy branch of a system with periodic 
boundary conditions in the x direction as analyzed in 
Ref. We identify this as a bulk mode. It corresponds 
to the first branch in Fig. El 

This is different than the situation we find when the 
boundary conditions given in Eq. (|10l) arc used. For com- 
parison, we show the square of the eigenfunction corre- 
sponding to the first branch in the upper panel of Fig. [6] 
for different values of the Poisson ratio a. We see that it 
is a constant function corresponding to a pure bulk mode 
only for a = 0, when the model correspond to a mem- 
brane. For finite values of a the eigenfunctions depend 
on the position over the ribbons and do not represent a 
bulk mode anymore. 




FIG. 6. Square of the normalized eigenfunctions f™(x) corre- 
sponding to the first (upper panel) second (middle panel) and 
third branch (lower panel) of Fig. [S] We also show for com- 
parison the square of eigenfunctions for ribbons with finite 
Poisson ratio a = 0.1 and a — 0.25. 

Otherwise the dispersion relation of the second branch 
could be well fitted by a law X(q) = q~\J clq + a\q 2 which 
disperses linearly near q = and then cuadratically for 
large values of q . The square of the eigenfunction cor- 
responding to these modes is shown in the middle panel 
of Fig. [S] We see that the most important distortion oc- 
curs when we approach the edges. This mode is therefore 
mostly an edge excitation and its behavior is qualitatively 
similar for both a = and a ^ 0. 

We arrive at the important conclusion: when a theory 
of a membrane is used, which wc assume is the correct one 




FIG. 7. (h {x,y)) obtained with the eigenfunctions corre- 
sponding to the second branch (upper panel) and third branch 
(lower panel) of Fig. [S] We also show, for comparison, the cor- 
responding quantities for a ribbon with finite Poisson ratios 
a = 0.1 and a = 0.25. 



for single- or few-layer graphene systems, the low energy 
phonon spectrum is decoupled into two branches with 
quite different physical interpretations. One corresponds 
to a bulk excitation and the other one is mainly an edge 
mode. 

Let us explore the consequences for the mean square 
amplitude of the height (h 2 (x,y)). Differently than in 
clamped ribbons, there are two gapless branches that give 
rise to a strong dependency with the length of the system 
~ L and a divergence when L — > oo. This can be inter- 
preted as an intrinsic instability of the system. However, 
it is known that the harmonic approximation is not valid 
for small values of q where the anharmonic coupling be- 
tween the in-plane and the out-of-plane mode a 10 ' 11 can 
stabilize the ribbon (or at least weaken the divergences) . 
Moreover, the external strains which are ubiquitous in 
real samples, will help to stabilize the system-^. 

The physical effect of the distortion on the ribbon will 
be quite different depending on the contribution of each 
of the two lower energy branches. The first one, with 
quadratic dispersion relation, is a bulk mode. We expect 
the same effect as the one analyzed in Ref. [l(| i.e. to 
produce an homogeneous rippling. More important for 
our goal in this paper is the effect of the second branch. 
In Fig. [7] we show the result for (h 2 (x,y)) calculated us- 
ing the eigenvectors and eigenvalues corresponding to the 
second (upper panel) and the third (lowest panel) branch 
of Fig. [5J As in the previous figure, we compare this re- 
sult for ribbons with finite Poisson ratios a = 0.1 and 
fj = 0.25 (using the values corresponding to the second 
and the third branches of its dispersion relation) . We see 
again that the maximum of (h 2 (x,y)} is at the bound- 
ary for each value of a. Regarding the contribution of 
the third branch, the correlation is strongly reduced as 
there is a gap for these excitations. A minimum around 
x ~ 0.28 is observed as well, corresponding to a node in 
the eigenfunctions. 
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FIG. 8. Correlation (ft(0.25, 0)/i(0.25, y)) obtained with the 
eigenfunctions and eigenvalues corresponding to the second 
branch (upper panel) and the third branch (lower panel) of 
Fig. [5] We also show, for comparison, the corresponding 
quantity for ribbons with finite Poisson ratios a = 0.1 and 
a = 0.25. The length of the ribbons is L — 1000 and its 
width W = 100. 



The height-height correlation function is shown in Fig. 
[5] as a function of y for fixed x = 0.25. Regarding the 
contribution of the third branch in the lower panel, we 
see a similar behavior than the one seen in the clamped 
ribbon analyzed in the previous section. The correlation 
function has very small values and it decays exponentially 
with the y distance. However, the contribution of the sec- 
ond branch is much more important and the correlation 
function decays much more slowly. This is a consequence 
of the absence of a gap for this mode. In fact, the ulti- 
mate value that this correlation assumes, increases with 
the length of the ribbon. 



Finally, let us briefly discuss the consequence of this 
result for a distortion of a free standing sample. It has 
been seen from early measurements that when graphene 
samples are put in a scaffold configuration, the free edges 
appear folded^. The fact that a scrolled configuration 
can be a stable deformation of a graphene edge has been 
previously shownii It is stabilized by an interplay be- 
tween the van der Waals and bending energies. The edge 
modes found in the present paper could be interpreted as 
precursive modes for this large distorted structure. That 
is to say, one possible sequence of events might be: first 
the edge oscillates rapidly with increasing amplitude and 
second it sticks to the rest of the sample due to the at- 
tractive van der Waals force. 



VI. CONCLUSION AND DISCUSSION 

In this paper we have studied the out-of-plane 
phononic spectrum and the height-height correlation 
functions of a graphene nanoribbon. Two different con- 
figurations were considered: clamped at the edges and 
free edges. When the ribbon is gripped, there are no true 
acoustic branches but a gap in the phononic spectrum. 
This gap leads to an exponential decay of the correlation 
and thermal excitations produce only local distortions of 
the crystalline order. This looks quite similar to the sit- 
uation in usual three dimensional solids. However, when 
the lateral dimension W increases enough, the character- 
istic correlation length also increases. This is coherent 
with the fact that in a square infinite membrane, the 
height-height correlation decays as a power lawii. 

When the edges are free, we find a quite remarkable 
decoupling between two different low energy phononic 
branches. Both of them go to zero energy for q —> 0. 
One disperses quadratically for all values of q. The eigen- 
vectors corresponding to this branch do not depend on 
the transversal coordinate of the ribbon. Therefore, in 
these modes the system is unaware of the existence of 
the edges. It is in fact a bulk excitation, similar to a 
lower excitation of an infinite membrane. The possible 
instability connected with this mode could lead the sys- 
tem to an homogeneous rippling. The other low energy 
branch disperses linearly for small values of q. The corre- 
sponding eigenvectors have their maximum values on the 
edges. It is a surface or edge mode. An instability is also 
associated with these modes. We claim that they could 
be the precursive modes resulting in folded edges as seen 
experimentally in suspended graphene samples. We re- 
mark that this decoupling between a bulk an edge mode 
is specific of a membrane-like model. In a model for a 
continuum plate in which a finite thickness is assumed, 
both modes are essentially edge modes. As graphene is a 
one-atom thick membrane, we expect that this fact would 
be specially significant for graphene nanorribons. 

A possible extension of our work could be the study 
of the interactions between the phononics modes and the 
conducting electrons in nanoribbons. It has been shown 
that fiexural phonons play an essential role as a limiting 
mechanism for the mobility in suspended graphene. We 
presume similar effects in nanoribbons. Moreover, when 
the edges of the ribbon are of zig-zag type, localized elec- 
tronic states near the edge are expected^. The phononic 
modes found in the present work should interact consid- 
erably with these confined electronic states. We expect 
them to play an important role in the transport proper- 
ties for these nanoribbons. 
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